Super-shell structure in harmonically trapped fermionic gases 
and its semiclassical interpretation 
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It was recently shown in self-consistent Hartree-Fock calculations that a harmonically trapped 
dilute gas of fermionic atoms with a repulsive two-body interaction exhibits a pronounced super- 
shell structure: the shell fillings due to the spherical harmonic trapping potential are modulated 
by a beat mode. This changes the "magic numbers" occurring between the beat nodes by half a 
period. The length and amplitude of the beating mode depends on the strength of the interaction. 
We give a qualitative interpretation of the beat structure in terms of a semiclassical trace formula 
that uniformly describes the symmetry breaking U(3) — > SO(3) in a 3D harmonic oscillator potential 
perturbed by an anharmonic term oc r 4 with arbitrary strength. We show that at low Fermi energies 
(or particle numbers), the beating gross-shell structure of this system is dominated solely by the 
two-fold degenerate circular and (diametrically) pendulating orbits. 

PACS numbers: 05.30.Fk 
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The realization of Bose-Einstein condensation in 
trapped dilute atom gases [l[ was a milestone in quan- 
tum physics, followed by a revolutionary development 
of both experiment and theory. Turning from bosonic 
to fermionic statistics, much current interest concerns 
the trapping and cooling of fermionic alkalis. The fo- 
cus presently often lies on pairing and the transition to 
a superfluid state ■ 

For these ultracold atomic fermi gases, not only is 
it possible to taylor the trap geometry, but also to ex- 
perimentally change the value of the scattering length 
for two-body collisions [f|. Being able to experimen- 
tally modify the interactions from attractive to repul- 
sive, entirely different interaction regimes can be probed: 
Atomic finite fermion systems are a unique laboratory to 
study fundamental quantum phenomena. 

It is apparent that these systems have much in common 
with atomic nuclei, with pronounced shell structures at 
low particle densities and a similar pairing mechanism for 
attractive short-range interactions U 0| ■ For weak repul- 
sive interactions between the trapped fermionic atoms, 
shell structure occurs in much analogy also to other fi- 
nite quantal systems, as for example electrons trapped 
in nanostructured semiconductor devices (so-called quan- 
tum dots |8j) or metallic clusters, in which the delocalized 
valence electrons are bound in the field of the metallic 
ions 0. 

Very recently we showed by self-consistent 

Hartree-Fock (HF) calculations that a harmonically 
trapped gas of fermionic atoms interacting by a weak 
repulsive two-body force may even exhibit a so-called 
super-shell structure. This means that the shell oscil- 
lations of the spherical harmonic oscillator are modu- 
lated by a beat structure, whereby the positions of the 
magic numbers are shifted by half a period between suc- 
cessive beats. Similar super-shell structure occurs in 
metallic clusters [llj]. Inspired by a semiclassical anal- 



ysis of Balian and Bloch in terms of the periodic or- 
bits in a spherical cavity [l2| , supershell beating patterns 
were predicted to occur in the abundance spectra of 
sodium clusters, and could later be observed experimen- 
tally 0. 

The spherical cavity model of Balian and Bloch [12| 
can, however, not be applied to the present system: For 
a dilute atomic fermi gas with short-range interactions 
trapped in a harmonic well, the semiclassical picture is 
different. We show that the super-shell structure origi- 
nates from the interference of diameter and circle orbits 
surviving the breaking of the U(3) symmetry of the har- 
monic oscillator by the leading anharmonicity term in 
the mean field. 

We confine the dilute gas of fermionic atoms by a 
spherical harmonic potential modeling an external trap 
@, interacting through a repulsive zero-range two-body 
potential ~ <z X}i<j " ( r -i ~ r j) where a is the s-wave 
scattering length. Due to the Pauli principle the S in- 
teraction only applies to fermions of pairwise opposite 
spin. We consider a fully unpolarized two-component 
system with two spin states, so that the total particle 
density is composed of two different densities of equal 
magnitude, n(r) = nJ(r) + n^(r) — 2rv(r). In the weak- 
interaction regime, the interaction energy density is given 
by grv (r)n-'-(r) = gn 2 (r)/4, where the coupling strength 
parameter g = Airft^a/m is introduced. This leads to the 
single-particle Hartree-Fock equation 



-— A + gn^(r) + V ho {r) 



V^(r)=^(r), (1) 



where Vh is the harmonic oscillator (HO) trap potential 
(for details, see Ref. [13]). In order to treat the interac- 
tion as a two-body process, diluteness of the gas requires 
that the interparticle spacing n -1 / 3 is much larger than 
the range of the interaction, and that na 3 <C 1. 

We solve Eq. ([T]) self-consistently on a grid under the 
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assumption of spherical SO (3) symmetry, which leads 
to states with (21 + l)-fold angular momentum degen- 
eracy. The HF interaction term is updated (with some 
weight factors) in each iteration according to gn'(r) = 

After convergence is obtained, the HF ground-state en- 
ergy of the A-particle system is summed up. In general, 
the ground-state energy as a function of A can be writ- 
ten as the sum of a smooth average part and an oscil- 
lating part, Etot = E av + E osc . The oscillating part, 
referred to as the shell-correction energy, or shell energy 
in short, reflects the quantized level spectrum {ei}. For a 
non-interacting Fermi gas in a spherically symmetric 3D 
harmonic trap, the leading-order term for the average en- 
ergy is found in the Thomas-Fermi approximation to be 
[j~6{ E>™ = (3N) 4 / 3 hw/4. For the repulsive interacting 
case, we find E av (g > 0) oc N a with a larger exponent 
a > 4/3. However, Eq. ([1]) with an interaction term lin- 
ear in the density is only valid for moderate g values and 
in practice we are close to a = 4/3 (e.g., a as 1.35 for 
g = 2). Contrary to the non-interacting case, and also to 
self-saturating fermion systems (such as nuclei and metal 
clusters) with a nearly constant particle density, it is not 
possible here to obtain the smooth part of the energy by 
a simple expansion in volume, surface and higher-order 
terms. We therefore perform a numerical averaging of 
the total HF energy over the particle number N in order 
to extract its oscillating part. 
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FIG. 1: The oscillating part of the ground state energy (in 
units of huj) as a function of 7V 1/3 for g = 0.2 (top, blue), 0.4 
(center, red) and 2 (bottom, green). The two lower curves are 
vertically shifted by —400 hu> and —600 hui, respectively. 

In the non-interacting case (5=0) the shell energy E osc 
oscillates with a frequency 27T3 1 / 3 as 9.06 as a function 
of N 1 / 3 and has a smoothly growing amplitude oc TV 2 / 3 . 
This follows from the exact trace formula [l6| for E osc of 
the 3D harmonic oscillator, whose leading-order term is 
given by 

~ (SAO^E cos (2nk(3N)^) . (2) 

fc=i 



Hereby k is the repetition number of the primitive clas- 
sical periodic orbit of the system with action Sq(E) = 
2ttE/lo. The argument of the cosine function in Eq. (jSJ) 
is simply k times Sq(E)/H, taken at the Thomas-Fermi 
value of the Fermi energy E F (N) = (ZNf^hw. The 
gross-shell structure is governed by the lowest harmonic 
with k — 1. 

Switching on the interaction, this scenario changes. A 
beating modulation of the rapid oscillations is found. In 
Fig. [T] we show the shell energy versus N 1 ! 3 for three 
values of the interaction strength, g=0.2, 0.4 and 2. A 
beating modulation of the amplitude of the shell energy, 
i.e., a super-shell structure, is clearly seen to appear for 
all cases. At small particle numbers and particularly for 
small g values, the shell energy is very close to that of the 
non-interacting system, given by Eq. ((2|). For larger inter- 
action strengths the super-shell structure is more clearly 
seen, and several beating nodes appear for g=2. With 
increasing interaction strength the amplitude of the shell 
energy oscillations becomes smaller. For example, for 
particle numbers around 80 3 as 500000, the amplitude of 
the shell energy is about 40 tko, which is only about 10 -6 
of the total ground-state energy. 

Through Fourier analysis of the calculated shell energy, 
two frequencies are seen to smoothly appear with increas- 
ing g value around the HO frequency (27T3 1 / 3 as 9.06). 
The super-shell features appear when the contribution 
to the effective potential from the interaction, gn\ is 
sufficiently large, i.e., at large values of g and A. We 
also observe that (almost) until the first super-node, i.e., 
A 1 / 3 as 28 in Fig. Q] (lower curve), the magic numbers 
agree with the HO ones (g — 0). Between the first two 
super-nodes, i.e., 28 < A 1 / 3 < 49 in Fig. [TJ the magic 
numbers for the interacting system are situated in the 
middle of two HO magic numbers, i.e., they appear at 
the maxima of the fast shell oscillations. Then, after the 
second super-node they roughly agree with the unper- 
turbed HO ones again. 

In the following we describe the major tools in an on- 
going semiclassical interpretation of these features [17| | . 
The U(3) symmetry of the unperturbed HO system is 
broken by the term SU = gn) in |T]), resulting in the 
SO (3) symmetry of the interacting system. The shortest 
periodic orbits in this system are the pendulating diam- 
eter orbits and the circular orbits with a radius corre- 
sponding to the minimum of the effective potential in- 
cluding the centrifugal term. These two orbits lead to 
the observed supershell beating [ijj • The above symme- 
try breaking had not been discussed in the semiclassical 
literature before. In a perturbative approach [l8j . it can 
be accounted for by a group average of the lowest-order 
action shift AS(o) brought about by the perturbation 
of the system: (ei AS ^°')oeU(3)- Hereby o is an element 
of the group U(3) characterizing a member of the un- 
perturbed HO orbit family (ellipses or circles). For the 
average it is sufficient to integrate over the 4-dimensional 
manifold CP 2 [19j, which for a perturbation 6U(r) = er A 
can be done analytically [13] • We therefore model the 
self-consistent numerical HF field by the following per- 
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turbed HO potential: 



valid for all strenghts e: 



V(r) = V ho (r) + -r\ 



(3) 



where the anharmonic term simulates the symmetry 
breaking effect of the part gn^(r) in ([T]). For small inter- 
action strenghts g, e is proportional to g. 

In the perturbative regime (s < 1) we have found the 
following perturbed trace formula for the level density: 



Sg pert (E) 



2eh 2 



-V 

7T ^- — ' 

k=l 



(-1)* 



kS ( 



kS d 



,(4) 



where kSd and kS c are the classical actions of the diam- 
eter and circle orbits, respectively. In the limit e — > 0, 



their difference goes as k(S c — Sd) 



kenEp 



(N)/co 5 



that ([3]) tends to the level density of the pure HO limit 
corresponding to J2])- 

To cover larger values of g and N, we have also devel- 
oped an analytical uniform trace formula [17f . Uniform 
here means that it reproduces the HO trace formula for 
the U(3) (higher symmetry) limit e = in a smooth way, 
analogously to those derived earlier for U(l) [2(| and U(2) 
symmetry breaking [2l| . In the large-energy (or large- N) 
limit, it yields the correct trace formula for the diameter 
and circular orbits, forming two-fold degenerate families, 



Sg(e) ~ i A i( e ) MkS d (e)/h) + A%(e) sin(kS c (e)/h)] . 
fc=i 

(5) 

Analytical expressions for the amplitudes and actions (in 
terms of elliptic integrals) are given in [17| • Eq. §5§ goes 
over into the perturbative trace formula (U in the limit 
e — > 0. This confirmes the statement made in that 
only the diameter and circle orbits are important in this 
limit. 

In the upper panel of Fig. [5J we compare the level den- 
sity for the potential ^ obtained with the perturbative 
trace formula ^ (dashed line, red) and the uniform trace 
fomula ([5]) (solid line, blue) for the value e = 0.005, us- 
ing only the lowest harmonics (k = 1). While there is 
good agreement up to the first super-shell maximum for 
this small value of e, for larger perturbations the two for- 
mulas only agree for very small energies. In the lower 
panel of Fig. [2] we compare the coarse-grained result of 
the uniform trace formula (solid line, blue) with the 
exact quantum-mechanical result (dashed line, red) for 
e = 0.01, demonstrating that the semiclassical interpre- 
tation of the super-shell beat in terms of circular and 
diameter orbits is perfect. 
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FIG. 2: Oscillating part Sg(E) of level density of the perturbed HO potential @ versus energy E (unit: hui), coarse-grained 
by Gaussian convolution over an energy range 7 = 0.5hu). Upper panel: comparison of perturbative trace formula @ (dashed 
line, red) and uniform trace formula ([5j) (solid line, blue) for the value e = 0.005. The dotted lines describe the amplitude 
for the unperturbed HO limit e = 0. Lower panel: comparison of uniform trace formula ([5} (solid line, blue) and exact 
quantum-mechanical results (dashed line, red) for e = 0.01. 
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At sufficiently high energies and perturbation 
strengths e, three-fold degenerate families of tori with 
rational ratios u> r : lu v — n : m > 7 : 3 of radial and 
angular frequency bifurcate from higher repetitions 
(k > 3) of the circle orbit [221 ] . The shortest of them, a 
star-like orbit with n : m = 7 : 3, is approximately 10 
times longer than the shortest diameter and circle orbits. 
Therefore, all the tori only contribute to finer quantum 
structures at higher energies. They can be included in 
the semiclasscial trace formula using standard techniques 

muni. 

The beat structure in E osc has some similarities with 
that found in nuclei @ and metal clusters Q • There are, 
however, two essential differences. 1. Those systems are 
self-saturating and have steep mean-field potentials that 
can be modeled by a spherical cavity [12|. The present 
system, in contrast, has a mean field with much smoother 
walls that are dominated at large distances by the confin- 
ing harmonic potential. 2. The super-shells in the cavity 
model come from the interference of the shortest periodic 
orbit families with three-fold degeneracy, as is usual in 
spherical systems 0, [24|. Here, however, the gross-shell 
structure comes uniquely from the diameter and circle 
orbits which are only two-fold degenerate, whereas the 
fully three-fold degenerate tori are so much longer that 
they only affect the finer quantum structures at higher 
energies. 



In conclusion, we have seen that the shell structure of 
weakly interacting fermions in a harmonic trap shows a 
pronounced beating pattern, with the single shell posi- 
tions changing by half a period length between the dif- 
ferent beat nodes. A Fourier analysis of the oscillating 
shell-correction part of the Hartree-Fock energy shows 
clear peaks at two slightly different frequencies. This we 
have interpreted semiclassically by the interference of the 
shortest periodic orbits generated by the breaking of the 
U(3) symmetry of the non-interacting HO system, which 
are the families of diameter and circle orbits. A more de- 
tailed quantitative interpretation using the uniform trace 
formula ([5]) derived in [17| is in progress. After extract- 
ing e and uj e ff from a polynomial fit to the numerical 
HF potential, we expect to describe the beat structure 
in the numerically obtained HF shell energies E osc (N) 
quantitatively in terms of classical periodic orbits. 
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